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Abstract. Using the newly developed VASP2WANNIER90 interface we have 
constructed maximally localized Wannier functions (MLWFs) for the Cg states of 
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the prototypical Jahn- Teller raagnetic perovskite LaMnOa at different levels of 
approximation for the exchange-correlation kernel. These include conventional density 
functional theory (DFT) with and without additional on-site Hubbard U term, hybrid- 
DFT, and partially self-consistent GW. By suitably mapping the MLWFs onto an 
effective Eg tight-binding (TB) Hamiltonian wc have computed a complete set of TB 
parameters which should serve as guidance for more elaborate treatments of correlation 
effects in effective Hamiltonian-based approaches. The method-dependent changes of 
the calculated TB parameters and their interplay with the electron-electron (el-el) 
interaction term are discussed and interpreted. We discuss two alternative model 
parameterizations: one in which the effects of the el-el interaction are implicitly 
incorporated in the otherwise "noninteracting" TB parameters, and a second where 
we include an explicit mean-field el-el interaction term in the TB Hamiltonian. Both 
models yield a set of tabulated TB parameters which provide the band dispersion in 
excellent agreement with the underlying ah initio and MLWF bands. 
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1. Introduction 



Perovskite transition-metal oxides challenge electronic structure theory since 
several decades, due to the variety of collective structural, electronic, and magnetic 
phenomena which are responsible for the formation of complex orbital- and spin- 
ordered states [H El [3]. A prototypical textbook example of this class of materials 
is the antiferromagnetic insulator LaMnOa. The ground state electronic structure 
of LaMnOs is characterized by the crystal-field induced breaking of the degeneracy 
of the Mn^"*" 3d^ manifold in the high-spin configuration {t2g)^{egY, with the t2g 
orbitals lying lower in energy than the two-fold degenerate Cg ones. Due to the strong 
Hund's rule coupling, the spins of the fully occupied majority t2g orbitals are aligned 
parallel with the spin of the singly occupied majority Cg states on the same site. The 
orbital degeneracy in the Cg channel is further lifted via cooperative Jahn- Teller (JT) 
distortions IHEl |6l |7] , manifested by long and short Mn-0 octahedral bonds alternating 
along the conventional orthorhombic basal plane, which are accompanied by GdFeOa- 
type (GFO) checkerboard tilting and rotations of the oxygen octahedra [H HI [10] (see 
figure [1]). As a result, the ideal cubic perovskite structure is strongly distorted into 
an orthorhombic structure with Phnm symmetry [HI |9], and a ci-type orbital-ordered 
(00) state emerges [11]. The corresponding occupied Cg orbital can be written as 
1^) = cosf|322_^2^^gi^||^2_^2^ [lailSl [Hills], with the sign of ^ ~ 108° alternating 
along X and y and repeating along z. This particular orbital ordering is responsible 
for the observed A- type antiferromagnetic arrangement below Tn = 140 K [161 IH]- It 
was found that long-range order disappears above 750 K, whereas a local JT distortion 
(without long-range order) remains (dynamically) active up to > 1150 K [6l [71 [Ti]. 

The question of whether the origin of orbital ordering should be attributed to a 
superexchange mechanism (0-mediated virtual hopping of electrons between nearest 
neighbor S = 2 Mn cations, associated with a local Coulomb electron-electron 
interaction: d'^d'^ ^ d\d^) [19] or to an electron-lattice coupling effect (structural- 
induced sphtting of the degenerate eg levels) [12] has been the subject of numerous 
studies [20l[2ll[22l|23l[2ll[25l[26l[ni[2Il[ll]. Considering that there is no clear 
experimental evidence to support one mechanism over the other, the employment of 
theoretical models and computer simulations has become an essential tool to explain 
the complicated coupling between structural and electronic degrees of freedom and 
to interpret the experimental observations. On the basis of model calculations, it 
has been recognized that the simultaneous inclusion of both superexchange and JT 
interactions is crucial to provide, to some extent, a satisfactory description of the 
observed transition temperatures Tn, Tqo and Tjt [201 [131 [H]. This approach typically 
relies on a suitable mapping between a realistic band structure calculated e.g. via 
density functional theory (DFT) [28] and an effective many-body Hamiltonian, which is 
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Figure 1. Representation of tlie JT/GFO distorted LaMnOa structure. Small red 
and light blue spheres indicate oxygen and manganese atoms, respectively, whereas 
the larger spheres refer to the La atoms. Plot done using the VESTA visualization 
program [iSj . 

often achieved by downfolding the relevant bands and constructing a locahzed Wannier 
basis [291 [131 EQl EH [32]. 

The quahty and characteristics of the Wannier representation inevitably depend on 
the underlying Kohn-Sham states. It is well known that the mean-field-type one-particle 
description of the electronic structure within the standard local density (LDA) [28] 
or generalized gradient (GGA) [33] approximations to DFT is incapable to correctly 
describe exchange and correlation effects in the so called strongly- correlated materials, 
resulting, among other failures [j| in much too small band gaps and magnetic moments [Ij. 
For this reason, the DFT-derived subset of orbitals is typically employed as reference 
for the one-electron (i.e. non-interacting) part of the effective Hamiltonian, where all 
approximated contributions coming from LDA/GGA exchange-correlation effects are 
subtracted in order to avoid double-counting p^. For example, in the DFT-I-DMFT 
method (combination of DFT and dynamical mean-field theory (DMFT)) [35], the 
effective Hamiltonian can be written as H = Hbft — H^c + Hjj, where -^dft is 
the Kohn-Sham Hamiltonian, i^^c accounts for the double-counting correction, and 
Hjj represents the Hubbard-like term which describes the electronic interactions in 
the strongly correlated bands. A critical issue of the DFT-I-DMFT approach is 
that a well defined expression for the double-counting potential is not known and 
several forms have been suggested [3H [36]. Karolak and coworkers have recently 
addressed this issue by treating the double-counting term as an adjustable parameter 

I We note that the underestimation of the band gap and related failures are of course also partly due 
to the intrinsic limitation of the Kohn-Sham approach, which is not meant to describe quasi-particle 
excitations correctly. 
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and suggested that the best agreement with experiment is achieved by setting the 
double-counting potential in the middle of the gap of the impurity spectral function |37] . 
Within this context, it is therefore justified to construct effective Hamiltonians starting 
from band structures obtained using different schemes, such as e.g. LDA+?7 [13] or 
hybrid functionals |38], which usually provide much better gaps for semiconducting 
materials than conventional DFT approximations and could therefore represent a more 
appropriate "non-interacting" reference for model calculations. 

For practical purposes, the most suitable starting point to study the physics of 
complex transition-metal oxides is probably the tight-binding (TB) scheme, which relies 
on a proper representation of the electronic structure using a localized basis set [H [3] . 
Some of the authors have recently shown that maximally localized Wannier functions 
(MLWFs) can be used to extract an effective TB description of the eg subspace in 
LaMnOs [3ll [32]. The calculated TB parameters can then be used to construct a 
simplified TB Hamiltonian in the form that is very often used for the description of 
manganites, ^tb = -^kin + -^Hund + -^jt + -^e-e, which then provides a very accurate 
representation of the underlying Kohn-Sham band structure. 

Motivated by the reasons outlined above, here we calculate MLWFs for LaMnOa 
using several different methods, including both the conventional GGA scheme and the 
more sophisticated GGA+U, hybrid functionals, and GW approaches. Besides providing 
a detailed description of the electronic and magnetic properties of LaMnOs at various 
levels of theory, we investigate how the corresponding differences in the treatment of 
exchange-correlation effects influence the specific features of the MLWFs and the TB 
parameters derived from them. 



2. Methodology and Computational Details 



2.1. DFT-based calculations 

All our calculations are based on DFT within the Perdew-Burke-Ernzerhof [33] (PBE) 
approximation to the exchange-correlation energy. The one-particle Kohn-Sham orbitals 
are computed within a plane-wave basis employing two different codes: (i) the program 
PWscf in combination with ultrasoft pseudopotentials included in the Quantum 
ESPRESSO package [39], and (ii) the projector augmented wave [IHIIIT] (PAW) based 
Vienna ab initio simulation package (VASP) [I2ll^. In particular, the PWscf program 
is used to benchmark the implementation of the VASP2WANNIER90 interface at PBE 
and PBE+f/ level. Due to the well known limitations of standard DFT in describing 
the electronic structure of "strongly-correlated" compounds, three different corrections 
to the PBE wavefunctions are adopted: (i) PBE+f/ : inclusion of a repulsive on-site 
Coulomb interaction U following the recipe of Dudarev et al. [B] ; (ii) Hybrid functionals: 
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suitable mixing between density functional and Hartree-Fock theory JIS] within the 
scheme proposed by Heyd, Scuseria, and Ernzerhof (HSE06, HSE hereafter) in which 
one quarter of the short-ranged exchange-correlation PBE functional is replaced by one 
quarter of the short-ranged part of Hartree-Fock exchange [IHl H?]; (iii) GW: explicit 
evaluation of the self-energy S = iGW within a partially self-consistent GWq procedure 
made up of self-consistent update of the eigenvalues in the Green's function G and a fixed 
screened exchange Wq, evaluated using PBE wavefunctions [HI HQ]. In accordance with 
previous studies [Ml Il9] , five iterations were sufficient to obtain quasiparticle energies 
converged to about 0.05 eV. 

These four methodologies (PBE, PBE+U, HSE and GW) differ in a few 
fundamental issues: (i) PBE relies on an approximate treatment of exchange- correlation 
effects; (ii) PBE+U contains the same PBE approximate correlation, but takes into 
account orbital dependence (applied to the d states of manganese) of the Coulomb and 
exchange interactions which is absent in the PBE; (iii) HSE includes a portion of non- 
local fully orbital dependent exact exchange and PBE correlation (iv) In GW exchange 
and correlation contributions are directly computed from the self-energy. 

An identical technical setup is adopted for VASP and PWscf calculations. All 
ground state electronic and magnetic properties are calculated for the experimental 
low temperature Pbnm structure reported in [8] using a regular F-centered 7x7x5 and 
6x6x6 k-point mesh in PWscf and VASP, respectively (reduced to 4x4x4 at the GWq 
level), and a plane wave energy cutoff of 35 Ry {^476 eV) and 300 eV in PWscf and 
VASP, respectively. Spin-polarized calculations were performed within a coUinear setup 
without the inclusion of spin-orbit effects. Except where otherwise noted, all PBE and 
PBE+f/ results discussed in the present work refer to PWscf calculations whereas HSE 
and GWq results are obtained using VASP. In both PWscf and VASP we include the 
Mn(3s), Mn(3p), La(5s), and La(5p) semi-core states in the valence. In PWscf the 
(unoccupied) La(4/) states are excluded from the ultrasoftpseudopotential, whereas 
they are present in the corresponding VASP PAW potentially 

To obtain the model TB parameters we perform additional calculations for a 
simplified crystal structure with the same unit cell volume as the experimental Pbnm 
structure, but which involves only the staggered (Q^-type) JT distortion and no 
GFO distortion and no orthorhombic deformation of the lattice vectors {Q^ = 0). 
See [501 ED E2] for more details and an exact definition of the different distortion 
modes. The amplitude of is 0.199 and 0.184 A in the experimental Pbnm and in the 
simplified JT{Q^) structure, respectively, and the amplitude of in the experimental 
Pbnm structure is -0.071 A . 

§ In the construction of the MLWFs within VASP we have shifted the La(4/) states to higher energies 
through the apphcation of a large U = 10 eV in order to avoid the overlap between La(4/) and 
unoccupied Mn(eg) states, which would otherwise deteriorate the disentanglement procedure. 
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2.2. Maximally localized Wannier functions 

A set of localized Wannier functions |w„t) corresponding to a group of bands 
that are described by delocalized Bloch states lipmk) is defined by the following 
transformation: 

N 

e-'^-^, (1) 



V f 
Iwut) = q / dfc 



.m=l 

where T is the lattice vector of the unit cell associated with the Wannier function, m is a 
band index, k is the wave- vector of the Bloch function, and the integration is performed 
over the first Brillouin zone (BZ) of the lattice. Different choices for the unitary matrices 
Uj''^ lead to different Wannier functions, which are thus not uniquely defined by ([T]). 
A unique set of maximally localized Wannier functions (MLWFs) can be generated by 
minimizing the total quadratic spread of the Wannier orbitals [5T] . 

Once the transformation matrices U}''^ are determined, a TB representation of the 
Hamiltonian in the basis of MLWFs is obtained: 

H= h^ZclT+^TCmT +h.c. , (2) 



TAT 



with 



nm 



(2i)! 



J BZ ; 



e-^'^ ' . (3) 



I 

Here, eik is the eigenvalue corresponding to Bloch function \ipik)- For cases where the 
bands of interest do not form an isolated set of bands but are entangled with other 
bands, a two step procedure for obtaining the unitary transformation matrices (which 
in this case are typically rectangular) is employed [52j. We note that T and AT in 
(IT])-(|3]) indicate lattice translations, whereas for crystal structures with more than one 
atom per unit cell, n and m generally represent combined orbital, spin, and site indeces, 
specifying the various orbitals at all sites within the primitive unit cell. 

Based on the projected densities of states (PDOS) calculated within DFT, we 
determine a suitable energy window for the construction of the MLWFs (more details 
follow in section 13. 2p . MLWFs are constructed by merging PWscf and VASP with 
the wannierQO code using the available PW2WANNIER90 tool [531 El] and the 
newly introduced VASP2WANNIER90 interface, respectively. Technical details on the 
construction of MLWFs within the PAW formalism can be found in Ref . [55] . Starting 
from an initial projection of the Bloch bands onto atomic Cg basis functions {Sz"^ — r^) 
and \x'^ — y^) centered at different Mn sites within the unit cell, we obtain a set 
of two Cg-like MLWFs per spin channel for each Mn site. The spread functional 
(both gauge-invariant and non-gauge-invariant parts) is considered to be converged if 
the corresponding fractional change between two successive iterations of the spread 
minimization is smaller than 10~^°. 
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Practical instructions for the use of VASP2WANNIER90: VASP uses wannier90 in 
library mode to generate all ingredients which are required to run the wannierQO code 
as a post-processing tool. 

Apart from the main wannierQO input file (wannierQO. win) the input files needed 
by wannierQO are [53]: (i) the overlaps between the cell periodic parts of the Bloch 
states (wannierQO. mmn), (ii) the projections of the Bloch states onto trial localized 
orbitals (wannierQO. amn), and (iii) the eigenvalues file (wannierQO. eig). This set of files 
is generated by VASP by setting LWANNIERQO = .TRUE, in the main VASP input file 
(INCAR). If the file wannierQO. win already exists, VASP will properly generate the files 
(i)-(iii) according to the instructions specified in wannierQO. win. If wannierQO. win does 
not exist, VASP will generate a default wannierQO. win file, which should be suitably 
modified in accordance to the keyword list described in the wannierQO user guide [56] 
in order to tell VASP what quantities to compute. Then, VASP has to be run again 
in order to create the additional wannierQO input files. To construct the UNK files 
(the periodic part of the Bloch states represented on a regular real space grid), which 
are required to plot the MLWFs, it is necessary to set LWRITE_UNK = .TRUE, in 
the INCAR file. In a spin-polarized calculation two sets of input files are generated 
(VASP2 WANNIERQO is employed only once to generate the files wannierQO. mmn, 
wannierQO. amn, and wannierQO. eig. These files are then used as input files for wannierQO, 
which is serially run for each energy window). Please refer to the online documentation 
of wannierQO for a detailed description of all relevant instructions [56]. 



3. Results and Discussion 



In this section we will first present and compare the electronic and magnetic 
ground state obtained within the various levels of approximation (PBE, PBE+U, HSE 
and GWq), before we will describe the downfolding of the resulting band structure 
by Wannier function decomposition. Finally, TB parameterizations corresponding to 
effective Cg models, either with or without explicit electron-electron interaction term, 
will be derived from these results, and implications of the different underlying band- 
structures will be discussed. 

3.1. Electronic and magnetic ground state 

The calculated band structures are displayed in figure |2] and the corresponding indirect 
{E{) and smallest direct (E^) band gaps are listed in table [TJ The calculated valence 
and conduction band spectra and the PDOS (corresponding to Mn(eg), Mn(t2g), and 
0(p) states), are represented in figure [3] and figure HI respectively. 

It can be seen from figure [2] that the eigenvalue dispersion in LaMnOa is 
characterized by an insulating state with indirect energy gap. By comparing with the 
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Figure 2. Calculated band structure along certain high-symmetry directions within 
the BZ. Each panel reports results obtained by a different method, as specified in the 
panel title. E=0 is aligned to the middle of the gap. 

Table 1. Collection of calculated (present work and previous studies) and 
experimental value for the indirect {E\) and direct (£'d) band gap of LaMnOa. 
The measured values refer to optical conductivity [63j [Ml [63 , Raman [66], and 
photoemission |67| experiments. 
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PDOS shown in figure HI it becomes clear that within all methods the Mott-Hubbard 
gap is opened between occupied and empty states with predominant Mn(ec,) character. 
While the width of the band gap differs strongly between the various methods, each one 
is in good agreement with previous LDA/GGA [53 EHl EHl [60] , (LDA/GGA)+[/ [Ml ED], 
and hybrid functionals [61], respectively (see table [1]). Our partially-self consistent GWq 
cannot be directly compared with the single-shot GqWo results of Nohara et a/. [62] 
since the latter depend much more on the initial LDA wavefunction and consequentially 
convey a smaller bandgap. 

Due to the inadequate treatment of exchange-correlation effects, conventional PBE- 
DFT leads to a significant underestimation of E^^^ = 0.75 eV compared to the 
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experimental values obtained from optical conductivity measurements (1.1 eV [53] , 
1.9 eV [61], 2.0 eV p]), Raman (2.0 eV p]), and photoemission data (1.7 eV [67]). 
In addition, the uppermost filled Mn(eg) bands (with energies in the region between 
—1.3 eV and 0.0 eV) are well separated from the lower-lying mostly Mn{t2g)- and 0{p)- 
like states (below —1.5 eV). In contrast, while the lower part of the group of bands 
immediately above the gap (up to about 2 eV) exhibits predominant local majority 
spin Cg character, these bands are strongly entangled with local minority spin states 
at slightly higher energies (between approximately 1-2 eV). The inclusion of the on- 
site interaction term within the PBE+U approach, separates these higher-lying local 
minority spin t2g states from the local majority eg bands directly above the gap for 
U > 2 eV. Furthermore, increasing U also increases the band gap (E™^+^ = 1.42 eV 
for U = 4 eV) and lowers the filled eg states relative to the bands with dominant Mn(t2c,) 
and 0(p) character, which leads to an appreciable overlap between these sets of bands 
around the F point for U = 4 eV. 

Changing to a more elaborate treatment of the exchange-correlation kernel, we 
observe that HSE provides a value of the bandgap {Ef^^ = 2.55 eV) that is significantly 
larger (by ~ 0.5 eV) than the experimental measurements. This is in line with 
previous hybrid functional estimates based on the B3LYP approach implemented within 
a Gaussian basis set [61]. By comparing the PBE and HSE band gap one could argue 
that a smaller portion of exact Hartree-Fock exchange should be included in the hybrid 
functional framework in order to obtain a better agreement with experiment. Indeed, 
a reduced mixing parameter Omix = 0.15 shrinks the direct gap down to 1.79 eV, 
almost on par with the photoemission measurements of Saitoh and coworkers [67], and 
with the more recent optical conductivity data of Jung et a/. [Ml EH], and Kriiger et 
al. [66] . LaMnOs therefore seems to represent another example for which the one-quarter 
compromise (mixing 1/4 of exact exchange with 3/4 of DFT exchange) is not the ideal 
choice [HE]. Finally, the parameter-free GWq technique leads to a quite satisfactory 
prediction of the band gap, E^'^'^ = 1.68 eV, and about significantly larger than the 
only previous single-shot (i.e. perturbative) GqWo study of Nohara et al. based on 
initial IDA wavefunctions [62]. Similarly to HSE and PBE+f/ (for U = 3 eV), GWq 
deliver eg bands around Ep well separated from the 0(p) and Mn(t2g) bands below 
and, to a lesser extent, above (there is an appreciable mixing of Mn(eg) and Mn(t2g) 
states along the T-Z-F path around 2 eV), in clear contrast with the PBE picture which 
predicts a certain degree of overlap between the eg bands and the higher lying t2g bands. 

In order to provide further assessment of the quality of the various methods 
in describing the electronic structure of LaMnOs we compare in figure [3] the 
simulated valence and conduction band spectra with the corresponding photoemission 
spectroscopy and X-ray absorption spectroscopy data [69] • For negative energies 
(occupied states) none of the four methods differs dramatically from the experimental 
spectrum, even though the multi-peak structures in the range of —7 eV to —4 eV seen 
within PBE+f/ and HSE do not have a clear experimental correspondence, whereas 
PBE and GWq profiles better follow the main three experimental peak/shoulders. The 
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Figure 3. Comparison between experimental [69] (blue squares) and calculated 
valence and conduction band spectra for PBE, PBE-|-t/ (C/ = 3 and 4 eV), HSE, 
and GWo . The calculated and measured spectra have been aligned by overlapping the 
valence band maxima and conduction band minima. 



situation is more critical for the unoccupied region, since none of the methods is capable 
to correctly reproduce the two-peaks structure characterizing the onset of the conduction 
band right above Ep. These two peaks could be interpreted as formed by Cg (lower one) 
and t2g (second ones) contributions and are described differently by the various schemes, 
following the corresponding band dispersions discussed in Fig. |2) (i) PBE both peaks 
merge in one single strong electronic signal, reflecting the large overlap between Cg 
and t2g bands right above Ep; (ii) in PBE+U the two peak are much too separated, 
reflecting the wide eg-t2g band splitting; (iii) HSE and GWq are rather similar. Their 
spectra are characterized by a lower Cg small bunch of states (onset of the conduction 
band spectra) associated to a more intense t2g-hke peak, but the GWq ^g/t2g splitting 
(si 1.4 eV) better matches the experimental one (~ 1.1 eV) as compared to the larger 
HSE splitting {k. 1.7 eV). From these results we can infer that GWq and HSE convey the 
most satisfactory picture in terms of peak position and corresponding spectral weight 
for both occupied and unoccupied states, with GWq better reproducing the splitting 
between the two lower conduction peaks. However, it should be noted that the relative 
weights of the two lower conduction peaks do not match with experiment, indicating 
that it is necessary go beyond the GW approximation to obtain a refined agreement 
with experiment, for instance using the Bethe-Salpeter equation (this is beyond the 
scope of the present study). We underline once more that unlike PBE+U and HSE 
(in which the proper adjustment of the parameters U and amix can cure the bandgap 
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problem and lead to values of the gap close to the experimental ones), the parameter- 
free GWq scheme is capable to provide a rather accurate picture without the need of 
any adjustable parameter. 

Next, we analyze the magnetic properties in terms of the nearest-neighbor magnetic 
exchange interactions within the orthorhombic ab plane (Jab) and along c ( Jc) [TOllGTlITT] . 
This will provide further insights into the performance of the various methods with 
respect to energetic properties of LaMnOs. By mapping the calculated total energies 
for different magnetic configurations onto a classical Heisenberg Hamiltonian H = 
— ^^ijLjJijSi-Sj, the following equations for Jab and Jc can be obtained (see also 

mm)- 

Efm — Eaaf = — 32 Jc (4) 
EcAF - Efm = 647^6 • (5) 

Here, -Efm corresponds to the total energy for the ferromagnetic (FM) configuration, 
whereas -Eaaf and Eqaf indicate the total energies associated with antiferromagnetic 
(AFM) ordering along z, and a two-dimensional checker-board like arrangement within 
the xy plane, respectively [6T| . 

The values of Jab and Jc obtained using the various methods considered within 
this work are listed in table [2] along with the calculated magnetic moments at the 
Mn site. We note that, due to the neglect of orbital degrees of freedom which in 
LaMnOs are strongly coupled to spin degrees of freedom, it is not obvious whether a 
classical Heisenberg model is well suited to give a complete picture of the magnetic 
properties of LaMnOs. Nevertheless it can still provide an accurate parameterization 
of the energy differences between the various magnetic configurations. However, the 
quantitative comparison with the experimental coupling constants derived from spin- 
wave spectra, i.e. small fluctuations around the AFM ground state, should be taken with 
care. In view of this, we can draw the following conclusions about the efficiency of the 
various DFT and beyond-DFT methods employed in the present study: (i) the magnetic 
energy differences exhibit appreciable variation between VASP and PWscf leading to 
differences of about 1-2 meV in the magnetic coupling constants. This is most likely due 
to the different pseudopotential technique employed in the two codes (PAW method vs. 
ultrasoft pseudopotential), which lead to qualitative differences especially at PBE+U 
level, as discussed below. A more elaborate discussion on the performance of different 
functionals and methods in predicting the magnetic couplings is given in Refs. [601 E] ) 
where it is concluded that the PAW values are very similar to the full potential FLAPW 
ones, (ii) In both codes PBE gives the correct A- AFM ground state, delivering a negative 
J^ (JJASP _ _2,i3 nieV, JfWscf _ _Q j^gY) ^ positive Jab {J^^^^ = 3.22 meV, 
jpwscf ^ 4_5Q j^gY)_ ^--^ rj.^^ uj^jjv correction to PBE decreases the Efm - Eaaf 
energy difference and eventually leads to the prediction of a FM ground state for U 
larger than a certain value. This critical value is rather different within the two codes 
used in this study: Jc becomes positive for ?7 = 2 eV and U = 4 eV, in PWscf and 
VASP, respectively. We note that this difference is almost entirely due to the difference 
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in the corresponding PBE results. The [/-induced changes in the magnetic couphng 
constant relative to the U = reference are nearly identical within the two codes, 
(iv) While Jc within HSE and PBE+f/ (VASP) are very similar for U between 2-3 eV, the 
ratio between and Jab is rather different within the two approaches, (v) Within the 
limitations regarding the applicability of a Heisenberg picture to LaMnOs stated above, 
HSE seems to be most consistent with the values of the magnetic coupling constants 
derived from neutron diffraction measurements of spin-wave spectra [72] and magnon 
data [75j. This further confirms the predictive power of HSE in describing exchange 
interactions in transition metal oxides, as compared to other available beyond-DFT 
schemes [73] . 

We can also see that all methods result in values for the local magnetic moments of 
the Mn cation that are within the range of variation of the experimental data. Generally, 
increasing U within PBE+U leads to a more localized magnetization density compared 
to PBE, and thus increases the local magnetic moments. 

On the basis of the above analysis both of the electronic and magnetic properties 
of LaMnOs, we can conclude that HSE and, when applicable, GWq (the calculation 
of magnetic energies at GW level to extract exchange coupling constants is presently 
not possible, or at least extremely difficult) are the most consistent with the available 
experimental data in terms of spectral properties, electronic structure and magnetic 
exchange interactions of LaMnOs. In view of this, we can now proceed to the 
discussion of the Wannier-based description of the Cg bands and the associated TB 
parameterization. 

3.2. Maximally localized Wannier functions 

In this section we present the details for the construction of the MLWFs with 
predominant Cg character from the calculated bands around the gap. In a TB 
picture, these MLWFs can be seen as "antibonding" bands resulting from the a-type 
hybridization between the Mn((i) and 0{p) atomic orbitals. Note that in this and the 
next section the discussion of the PBE+U results refers to the representative value of 
U = 3 eV, unless explicitly stated otherwise. 

Figure S shows the PBE and beyond-PBE (PBE+U, HSE and GWq) band 
structures and the corresponding PDOS with Mn(eg), Mn(t2g), and 0(p) character. 
Apart from the obvious hybridization between Mia{d) and 0{p) states, "Cg-like" orbitals 
at a certain site can hybridize with "t2g-like" orbitals at a neighboring site as a result 
of the tilt and rotation of the oxygen octahedra. This leads to bands with mixed 
^g/hg character (note the bands around the gap with strong PDOS components of both 
Cg and t2g character). Due to this strong mixing it is not possible to construct 8 e„ 
character MLWFs within one energy window used in the disentanglement procedurel]]! 

II In the antiferromagnetic case each band is of course two- fold degenerate with respect to the global 
spin projection. Here and in the following we refer to such pairs of spin-degenerate bands as "one 
band" . 
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Table 2. PBE, PBE+C/, HSE and GWq derived magnetic exchange parameters (meV) 
and magnetic moment at Mn sites fi (pb)- The experimental and previously published 
computed data are taken from: Ref. [60], Ref. |6l], " Ref. [72], Ref. [8], " Ref. [74], 
and f Ref. 75, 



PWscf 



PBE 


4.56 


-0.81 


6.01 


U = 2 eV 


5.02 


0.37 


3.82 


U = 3eY 


5.30 


0.98 


3.89 


[/ = 4 eV 


5.63 


1.55 


3.96 




VASP 




PBE 


3.22 


-2.13 


3.50 


[/ = 2 eV 


3.54 


-0.84 


3.68 


[/ = 3 eV 


3.57 


-0.30 


3.76 


[/ = 4 eV 


3.61 


0.17 


3.83 


HSE 


2.56 


-0.53 


3.74 


GWq 






3.51 




Previous studies 




GGA+i7 ([/ = 2 eV)" 




-1.30 


3.46 


B3LYP'' 


2.09 


-1.01 


3.80 


Expt 




-1.16^ 


3.87^ 3.7±0, 




1.67^ 


-1.21^ 





The corresponding energy window would inevitably also contain the local minority 
spin '%g" bands. Since due to the GFO distortion these bands can hybridize with 
the minority spin "cc," bands, this would lead to MLWFs with strongly mixed eg/t2g 
character. To circumvent this problem, we therefore construct two separate sets of 4 
local majority and 4 local minority spin MLWFs using two different energy windows |3T] . 
These energy windows have to be chosen carefully for each individual method. (This 
problem is not present for the purely JT(Q^) distorted structure, from which we derive 
most of the model parameters, see section 13. 31 In this case we calculate a full set of 8 
MLWFs). 

To find a suitable energy window is quite straightforward for the local majority 
spin case. The upper bound of the energy window is determined by the upper bound 
of the highest (in energy) peak of the local majority spin Mn{eg) PDOS, while the 
lower bound of the energy window should be placed above the occupied bands with 
strong 0{p) and/or local majority spin Mn(t2g) character. It can be seen from figure H] 
(and perhaps more clearly from figure [2l), that both the lower and the upper bound 
fall within small gaps separating the bands within the energy window from other bands 
at lower and higher energies. Furthermore, for PBE+U and HSE the MLWFs can be 
constructed from a completely isolated set of bands, whereas in the case of PBE and 
GWq additional bands with predominant minority spin Mn(t2g) character are included 
in the energy window. However, due to the different local spin projection, these latter 
bands have no noticeable effect on the final MLWFs. 
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Figure 4. Effective Cg MLWF bands (thick red lines) for LaMnOs superimposed to 
the ah initio electronic bands (gray thin solid/dotted lines) and associated normalized 
PDOS (to the left and right of the band structure plots) corresponding to Mn(eg) (red 
filled areas), Mn(t2g) (green lines), and 0(p) character (blue dots). In the left /right 
PDOS graphs, Mn((i) PDOSs correspond to the local majority/minority Mn sites while 
the 0(p) PDOS is calculated as an average over all O sites. The two energy windows 
used in the wannier-downfolding are indicated by dashed and dot-dashed lines. The 
Fermi level (E=0 eV) is set in the middle of the gap. 



For the local minority spin MLWFs, the upper bound of the energy window can be 
found in the same way as for the local majority spin bands. Within PBE the lower bound 
is also easily determined, since it falls within a small gap separating the local minority 
spin bands with predominant Cg and character. However, no such gap exists within 
PBE+f/, HSE, and GWq, and it is thus not possible to fully exclude the t2g character 
from the resulting MLWFs. Instead, the lower bound of the energy window has to be 
carefully adjusted by manually checking the Cg character of the calculated MLWFs in 
real space. 

The band dispersion of the so-obtained MLWFs is shown in figure H] as thick red 
lines. The 4 (energetically lower) local majority MLWF bands follow very closely the 
underlying PWscf/VASP bands and the overall dispersion is very similar for all methods. 
Despite the strong band-entanglement, the dispersion of the 4 (energetically higher) local 
minority MLWF bands is also very similar within all methods. Only the energetically 
lowest local minority spin band within PBE+f/ and HSE exhibits strong deviations from 
the corresponding PBE and GWq case. This is due to the above-mentioned difficulty 
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to exclude the t2g character in a controlled way. Conclusions drawn from such sets of 
MLWFs should therefore be taken with care. Overall, we note that the similarities in 
the band structure and PDOS between PBE and GWq as well as between PBE+U and 
HSE, regarding the degree of hybridization between Mn(eg), Mn(t2g) and 0{p) orbitals, 
that have been pointed out in the previous section, are also reflected in the MLWF 
bands. 

To further demonstrate the similarities between MLWFs calculated at different 
levels of theory, we show in figure |5] the real space representation of the 2 MLWFs 
localized at a certain Mn site, projected on the xy plane. The dominant Cg character 
at the central Mn site together with the "hybridization tails" of mostly p character at 
the surrounding O sites is clearly visible for all MLWFs and methods. For the local 
majority spin MLWFs (1st and 3rd row), there is essentially no visible difference in 
orbital character between PBE and PBE+[/, only the 0{p) tails are marginally stronger 
if the Hubbard U correction is applied. At the HSE level, both local majority MLWFs 
exhibit significant x/y asymmetry, leading to more pronounced 0{p) hybridization tails 
along the short and long Mn-0 bond for the |3z^ — r^)-like and |a;^ — ?/^)-like function, 
respectively. Within GWq, the central e^-like part as well as the 0{p) tails are less 
asymmetric than for HSE, and appear similar to PBE/PBE+?/ for both local majority 
MLWFs. In comparison with the local majority MLWFs, the 0{p) hybridization tails 
of the local minority MLWFs (2nd and 4th row) are generally less pronounced. There 
is no significant difference between the local minority spin MLWFs calculated using the 
different methods. Even at the PBE+U and HSE levels, for which the admixture of the 
t2g character could not be controlled systematically, there is no apparent difference in 
comparison with PBE and GWq. The orbitally ordered states resulting from this set of 
MLWFs basis set is shown in Figj6] in terms of charge density isosurfaces of the highest 
occupied and lower unoccupied orbitals associated to the Cg bands below and above Ep 
in the lower energy window as defined in Fig. HI This plot clearly show the staggered 
ordering at neighbouring Mn sites and the significant p — d hybridization at the oxygen 
sites. As a comparison we provide in Fig. [7] the corresponding staggered ordering 
associated to the highest occupied e^-like bands as obtained from the full ab initio self- 
consistent charge density (without downfolding) within the various methods employed 
in the present study. The similarities between the ab initio and wannierized orbital 
ordering is a further demonstration of the quality and reliability of our wannierization 
procedure. 

3.3. Tight binding model Hamiltonian 

As already mentioned in the introduction, the electronic Hamiltonian of the eg manifold 
in manganites is generally described within the TB formalism as a sum of the kinetic 
energy H\an and several local interaction terms, the Hund's rule coupling to the 
t2g core spin -f^Hund, the JT coupling to the oxygen octahedra distortion -ffjT, and 
eventually the electron-electron interaction H(,-c^ which can be written as (see e.g. 
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Figure 5. Real space representation of the four Cg MLWFs corresponding to a certain 
Mn site, projected on the xy plane cutting through the Mn site. Black iso-lines 
correspond to ±N/^/V with integer > 1, the white region is defined by values in 
the interval [~l/^/V , + l/\/V], where V is the volume of the unit cell. Blueish/reddish 
hue denotes negative/positive values of MLWFs and Mn and O atoms are shown as 
blue and red spheres, respectively. 



PRE 



PBE+ (7 



HSE 



CW, 



(a) 



(b) 



mioccupiod 

V ^ 

V 



* * ** * 

^. 

r or J 


ft* ♦ . %• *• 


%* %• •# 




^ mi 


* > • 


fir 

• 

« « * • 


^ % ' % 

/■ /■ 

^ # % # 

^ % * % 


• » t Ik 



Figure 6. Charge density isosurfaces of the orbitally ordered states associated to the 
highest occupied (a) and lower unoccupied (b) MLWFs orbitals. Color coding and 
symbols are the same as in Fig. [S] 
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Figure 7. Charge density isosurfaces of the highest occupied Cg orbitals (from Ep 
to the lower energy bound as defined in FigH]) showing the orbitally ordered state of 
LaMnOa obtained using the different methodologies employed in this study. 



Refs. [31 [231 1501 [271 EB Eg) 

Hkm = - 5Z 4,a{ii+AH)^'^,«{J?+AK)6{R)Ca,6(H) , (6a) 
a,b,R,AR,a 

H-Rund = —Ju ^ Sr ^ C^,a(H)'''o''^'^o''.«{-R) ' (6^) 
R a,(j,a' 

a,b,R,i,a 

^e-e = 1 ^'^''^dcl^aiRfi' ,b{Rf'y' AR)Ca,c{R) ■ (6cO 

a,b,c,d,a,a' 

Here, c^^a{R) and the annihilation and creation operators associated with 

orbital \a) and spin a, centered at site R. Furthermore, t„^a{R+AR)b{R) are the hopping 
amplitudes between orbitals at site R and R + AR, r^^ are the standard Pauli matrices, 
Jh is the Hund's rule coupling strength, Sr is the normalized core spin at site 
R, X is the JT coupling constant, and is the amplitude of a particular JT mode 
{i = {x,z}). In our TB analysis we will only consider the electron-electron interaction 
within a mean-field approximation and use a simplified version of Eq. ( l6dl) corresponding 
to Uaaaa = Uabab = ^Av and all other interaction matrix elements set to zero, which is 
consistent with the PBE+f/ treatment according to Dudarev et al. [44]. The resulting 
shift in the one-electron potential due to the electron-electron interaction then becomes 

Va,ab = i^^ab — fT'u,ab) , (7) 

where U-w is the Hubbard parameter in the basis of MLWFs and Ucr^ab are the 
corresponding occupation matrix elementsl^l 

The model parameters (t„^a{R+AR)b{R), Ju, A, Uw) which determine the TB model 
Hamiltonian can in principle be obtained from the Hamiltonian matrix elements in 

% Here and in the following we often suppress either site or spin indeces or both of them, unless the 
corresponding values do not become clear from the context. Apart from the hopping amplitudes all 
quantities are diagonal in site index. In addition, for the collinear configurations of core-spins Sr 
considered here, the Hamiltonian and all quantities involved are also diagonal with respect to the 
global spin projection. 
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the MLWF basis. We note that AT in ([2]) refers to lattice translations whereas AR in 
(16 gp refers to the relative position with respect to the lattice of Mn sites. In the following 
we will therefore use the following simplified notation: — )■ h^^^^j^, — )■ where 
Ai2 = R' — R + AT. Then a and b correspond to the two effective Cg orbitals centered 
at individual Mn sites separated by AR. In order to further simplify the notation for 
the hopping amplitudes, we choose one Mn site as the origin [R = 0) and align the x 
and y axes of our coordinate system with the directions corresponding to the long and 
short Mn-0 bond of the JT(Q^) mode, respectively. We then define the vectors x, y, z 
according to the nearest-neighbor spacing of the Mn sites along the respective axes. 

Our TB parameterization is based on the procedure described by some of the 
authors in [31], with certain modifications, explained in the following. In |;31j it was 
shown that, at least at the PBE level, the influence of an individual structural distortion 
(JT or GFO) on the Hamiltonian matrix elements /if^^ expressed in the basis of e^-like 
MLWFs is to a great extent independent from the other distortion, and that furthermore 
the magnetic configuration has only a weak influence on the resulting model parameters. 
The TB parameterization was therefore based on various model structures with both 
FM (which always leads to a metallic system) and A-AFM order, with individual 
structural distortion modes frozen in. Due to the significantly increased computational 
cost of the HSE and GWq methods in comparison with PBE (in particular for the 
metallic state for which a dense k-points mesh is required to achieve a well converged 
solution), it is desirable to derive the TB parameters from as few (and if possible 
insulating) model structures as possible. In the present study, we therefore construct 
the TB parameterization from only two crystal structures: the purely JT(Q^) distorted 
structure and the experimental Pbnm structure, in both cases with A-AFM order, which 
then yields an insulating solution. As we will show in table |3l the TB parameters derived 
in this way at the PBE level deviate only marginally from the parameters found in the 
previous study [31] . 

In the following we describe the modified method we use to construct parameters 
of the TB model ( l6aj) -( l6dl) . Many of the simplifications on which our effective TB 
description of LaMnOs is based can be understood from the MLWF matrix elements 
shown in figure [HI and will be discussed in the remainder of this section. We will first 
consider an effectively "noninteracting" case in which we neglect the term i?e-e and 
consider how the more sophisticated beyond-PBE treatment of the exchange-correlation 
kernel affects the hopping, JT, and GFO-related parameters. We name this approach 
Model 1. Then, we discuss an alternative way which involves an explicit treatment of 
He-e in the model Hamiltonian within mean-field approximation. This allows us to 
obtain estimates for the corresponding on-site interaction parameters, by keeping the 
conventional PBE description as reference. We call this Model 2. Further technical 
details can be found in the Appendix. 

3.3.1. TB parameterization with implicit el-el interaction: Model 1. As shown in [31] 
for the PBE case, good agreement between an effective Cg TB model and the underlying 
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Kohn-Sham band structure can be achieved by considering hopping only between nearest 
neighbor Mn sites, next-nearest Mn neighbors, and second-nearest Mn neighbors along 
the X, y, and z axes, described by parameters t^*, t^^, and t^^, respectively (see 
Appendix). Thereby it is necessary to take into account the spin dependence of the 
nearest neighbor hopping amplitudes. This can also be seen from figure [H]( a) and (b), 
where (for PBE) the difference between (/i^a)^ ^^"^ i^aa)^ (from which t'^'^ and are 
calculated using (11. 2p ) is indeed significant. On the other hand, the further neighbor 
hoppings {t^y and t^^) show only negligible spin dependence, and are therefore calculated 
from the corresponding spin averaged Hamiltonian matrix elements. We note that s in 
t^'^ should be read as a local spin index (i.e. relative to the orientation of the local 
core-spin Sr) corresponding to the sites between the electron hops, which can have the 
values ±1 corresponding to t/i- The parameters and thus represent hopping 
amplitudes within FM ordered planes. As a result of the GFO distortion, and are 
reduced by a factor {1 — rj^), where rj^ is determined from the ratio of the t^'' calculated 
for the Pbnm and JT(Q^) structures (see (11. 3p in Appendix). The hopping amplitude 
t'^'^ between A-AFM ordered planes is then calculated as average of and t^'^. As 
also shown in [31] , the JT distortion induces a strong splitting between the nondiagonal 
elements of the nearest-neighbor hopping matrix within the xy plane (see the differences 
between /i^2 ^21 figure [8](a,b)), which is parameterized via a non-local JT coupling 
strength A (see (11. 5p in Appendix). 

Within model 1 only two contributions to the on-site part of the TB Hamiltonian 
are considered: the Hund's rule coupling i^Hund and the Jahn- Teller coupling Hjt. The 
strength of the Hund's rule coupling Jh is determined from the spin splitting of the on- 
site diagonal matrix elements h^a for the Pbnm structure, averaged over both orbitals 
(see Eq. (II. 8p ). The JT coupling strength A'^ for local spin-projection s is determined 
according to Eq. fl6cP from the splitting of the eigenvalues of the on-site Hamiltonian 
matrix and the JT amplitude for the purely 3T{Q^) distorted structure. As can be 
seen from figure [8](c,d), the corresponding matrix elements are strongly spin-dependent, 
leading to large differences in the corresponding JT coupling constants. Similar to the 
hopping amplitudes. A** is reduced by a factor (1 — rjl) due to the GFO distortion, which 
is determined from the ratio between A* calculated for the Pbnm and JT{Q^) structures. 

Table [3] lists the obtained TB parameters corresponding to Model 1 calculated 
within the various levels of approximation. Both hopping amplitudes and JT coupling 
strength correspond to the case without GFO distortion. It can be seen from the 
first two rows of table [3] that the parameterization we use in the present study yields 
only marginal differences for the PBE hopping parameters and Hund's rule coupling 
in comparison with [31]. This corroborates the quality of our TB parameterization 
based on only two structures (JT(Q^) and Pbnm with A-AFM order). Note, that here 
we use a crystal structure derived from low-temperature measurements [8], whereas in 
[31] the room temperature measurements of Ref. p] have been used. The JT coupling 
parameters differ slightly more from [31], due to the revised definition of A** used in 
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Figure 8. Hamiltonian matrix elements in the basis of MLWFs for the experimental 
Pbnm structure: nearest-neighbor terms corresponding to local majority (a) and 
minority (b) spin projection, diagonal (c) and off-diagonal (d) on-site terms. Local 
majority and minority spin projections are indicated by up and down triangles, 
respectively. Left /right parts of the horizontal axis corresponds to PWscf/VASP 
results. 



Table 3. The TB model parameters as derived from PBE and beyond-PBE band 
structures (Model 1, in PBE-I-U we used U=3 eV). Since the PBE-I-C7 values of r]f 
and are unreliable (see text), we use the corresponding PBE values (in brackets) 
to compute the TB bands displayed in figured Units: t^^, t^^, t^^ , t"^^ in meV; A in 
meV/A; Jh in eV; A''", in eV/A; rjj , r]f, rjx are unit-less. 
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HSE 


750 
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50 


0.40 


0.20 


2.42 


10.25 


0.96 


0.28 


GWo 


746 


469 


490 


13 


50 


0.24 


0.41 


1.90 


4.43 


0.88 


0.04 



the present study. Another important change arises from the use of 3 separate GFO 
reduction factors rji, ri^, and ri\, instead of using one averaged value as it was done in 
[31]), which provides a more accurate TB description of the MLWF bands. It can be 
also seen from table [3] that at the PBE level, there is essentially no difference between 
the hopping amplitudes calculated using PWscf and VASP. There is a 12 % difference 
in Jh between PBE(VASP) and PBE(PWscf), which could be related to the noticeable 
differences in the energetics of the various magnetic configurations discussed earlier. 
Comparing the parameters obtained from the beyond-PBE methods with the pure 
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PBE case, we observe that the hopping parameter is generally increased in all beyond- 
PBE methods. As was shown in [32], this can be understood within an extended nearest 
neighbor TB model including both Mn{d) and 0{p) states, from which an effective Cg- 
only model can be derived in the limit of large energy separation Sdp between the d and 
p orbitals. The effective hopping in the Cg model is then given in terms of the nearest 
neighbor hopping amplitude t^p of the extended d-p model as t^J^ = t'^p/£dp- The increase 
of t^^ is therefore consistent with the observation that all beyond-PBE methods lower 
the Cg bands relative to the lower-lying oxygen p bands. The small decrease of t'^^ within 
PBE+f/ (for small values of f/ < 2 eV) can be explained in the same way, since here the 
corresponding energy separation between 0(p) and Mn(eg) increases. The JT parameter 
A is generally very similar for PBE, PBE+f/, and GWq, while a strong enhancement 
of A can be seen for HSE, which is consistent with the strong x/y asymmetry of the 
corresponding MLWFs seen in figure |5]^c). Since the changes of the already rather small 
further-neighbor hoppings within the beyond-PBE methods are very small, we use the 
corresponding PBE values for simplicity. The GEO reduction factors for the hopping 
amplitudes, 1]} and rjj, are slightly decreased within GWq, whereas rjj is increased for 
PBE+f/ and HSE, and rjj is strongly decreased in HSE. Due to the strong mixing 
between minority spin Cg and t2g bands within PBE+f/, which was already discussed 
in section [3l2] (see also figure |9]^b)), the determination of rij is rather unreliable in this 
case, and we therefore use the corresponding PBE value. We note that the same effect 
also leads to the strong changes in the 1 cal minority hopping matri xelements within 
the xy plane calculated within PBE+f/ for f/ > 3 eV (see figure [Hl^b)). Using the HSE 
and GWq methods we do not encounter this problem. 

For all beyond-PBE methods, a significant increase of Jh and A^ can be observed, 
which in the TB model gives rise to an increase of the spin splitting and the band 
gap, respectively. The change of A"^ compared to PBE is very small for both HSE and 
GWq. Due to the inaccurate treatment of the minority spin bands, PBE+f/ gives an 
unrealistically small value of A^*" = 0.30 eV/A, which we therefore substitute with the 
corresponding PBE value. While r]\ does not change significantly for small values of the 
Hubbard U, a small increase (significant decrease) is observed for HSE (GWq). 

To assess the quality of our parameterization we now use the TB parameters 
tabulated in table [3] to compute the resulting Cg band structure. In figure [9]^a) and (c), 
we compare the band dispersions of the TB model (blue filled circles) and the MLWFs 
(thick red lines) for the experimental Phnm structure within the PBE approximation. 
Despite the many simplifications made in the construction of the model parameters, the 
TB model can reproduce the MLWF bands to a remarkable accuracy (for both PWscf 
and VASP). The reliability of the beyond-PBE TB representation can be appreciated 
by the overall excellent match between the TB and MLWFs bands shown in figure |9]^b), 
(d) and (e), which exhibit the same quality as observed at the PBE level. This is 
particularly true for the band gap, whose method-dependent changes (see table [1]) are 
perfectly reflected in the TB descriptionE] 

+ The MLWF and TB bands were aligned by minimizing a mean deviation which was calculated as 
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Figure 9. Comparison of the band dispersion corresponding to MLWFs (red lines), 
the TB Model 1 using parameters given in table |3] (blue circles) , and the TB Model 2 
with interaction parameters given in table |4] (green crosses) . 

3.3.2. TB parameterization with explicit el- el interaction: Model 2. Now, we turn 
our attention on the alternative TB parameterization in which we attempt to treat 
the modifications induced by the beyond-PBE methods as a perturbation to the 
"noninteracting" PBE description by exphcitly considering the el-el interaction ( l6dj) 
and using the simplified mean-field approximation ([7]) in the TB Hamiltonian. It is 
clear from the discussion in the preceeding section that it is not straightforward to 
parameterize the hopping amplitudes in terms of t/w. We will therefore limit ourselves 
to analyzing the effect of ([7]) on the local Hamiltonian, which can be represented as 2x2 
matrix in terms of the two local Cg states in the following form: 

£!ocai = lo - Uy,n^ (8) 

with 

= 1 {\Uy^ -Jh-s)- yQ\^ - X'QY ■ (9) 

By identifying (jH]) with the corresponding MLWF matrix, we obtain the local spin 
splitting as a combination of Hund's rule coupling and el-el interaction: 

W - {hlY = U^\nl - nl) + 2Jr"^ , (10) 

from which we can calculate t/-^'' by averaging over the two orbital characters and 
using the previously determined PBE value for the Hund's rule coupling. Thereby, the 
occupation matrix elements in the basis of MLWFs are calculated as 

nnm=rde[ dkJ2(u^::A*6ie-e,,)u!^\ (11) 

J-oo JbZ I ^ ' 

where is the Fermi energy. 

an average of the corresponding eigenvalue differences over all bands and k-points. The maximum and 
mean deviation is very similar for all methods and does not exceed 0.37 and 0.12 eV, respectively. 
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Table 4. The interaction parameters determined in Model 2. Note, that in Model 
2 the on-site parameters are set to the PBE values while the hopping parameters are 
set to the values given in table [3] Units: all quantities are in eV except An''' which is 
unit-less. 



Interaction parameters 























PWscf 








PBE 


1.56 


1.09 


0.71 








PBE+C/ 


2.16 


1.66 


0.80 


2.40 


0.70 


0.42 








VASP 








PBE 


1.33 


1.04 


0.70 








HSE 


2.42 


3.10 


0.89 


4.37 


2.31 


0.51 


GWo 


1.90 


1.80 


0.70 


2.30 


1.09 


0.30 



In a similar way we can obtain another estimate for the Hubbard parameter, 
from the total JT induced splitting within the majority spin eg orbital manifold, 
expressed through the difference in eigenvalues of the local Hamiltonian: 

Ae^ = 2A^(P^^)V(Q^)' + (Q^)' + t^w = Ae^^^^^) + f/^^An^. (12) 

Here, An^ is the difference in majority spin eigenvalues of the MLWF occupation matrix 
and we have used the observation that, to a very good approximation, both and n 
can be diagonalized by the same unitary transformation. The difference between the 
corresponding transformation angles is less than 0.6° for the Pbnm structure. Since 
the difference is somewhat larger for the JT(Q^) structure (up to ~ 6°) we derive the 
interaction parameter from the MLWF Hamiltonian of the Pbnm structure. The 
resulting values of and are given in table HI 

It can be seen that within PBE+U, the parameter is almost as large as the 
value of f/ = 3 eV used for the Hubbard parameter within the PBE-|-?7 calculation, 
whereas the parameter is significantly smaller than that. We note that, as discussed 
in [32], the Hubbard correction within PBE-|-?7 is applied to rather localized atomic- 
like orbitals, whereas the parameter f/w corresponds to more extended e^-like Wannier 
orbitals. The JT splitting is strongly affected by hybridization with the surrounding 
oxygen ligands and is thus quite different for atomic-hke and extended Wannier states 
[32] . As a result, is quite different from the U value used within FBE+U, and the 
smaller value of can thus be related to the fact that the electron-electron interaction 
is more screened in the more extended effective Cg Wannier orbitals. On the other hand, 
the similarity between and the U value used within FBE+U indicates that the 
local spin-splitting is more or less the same for both sets of orbitals, which is consistent 
with the view that this splitting is essentially an atomic property. A similar difference 
between and is also observed for HSE and GWq. The large values of f/w 
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delivered by HSE reflects the larger spin splitting and band gap in the corresponding 
band structure compared to PBE+f/ and GWq. 

The large difference between the two parameters U^"^ and U^^^ also indicates that 
it is not possible to map the electron-electron interaction effects manifested in the 
on-site matrix corresponding to effective Cg orbitals to only one interaction parameter 
while using PBE as "noninteracting" reference. Similar conclusions have already been 
reached in [32] for the PBE+f/ case. From the current study we can conclude that 
the modification of the local spin splitting (described by U^'^) and the enhancement 
of the JT induced orbital splitting (described by U^^) that arise in the Kohn-Sham or 
GWo quasiparticle band structures due to the beyond-PBE treatment of exchange and 
correlation, are not compatible with a simple mean-field Hubbard-like correction to an 
otherwise "non-interacting" TB Hamiltonian with two effective eg orbitals per Mn site 
and only one parameter describing the electron-electron interaction. This leads to an 
important conclusion of the present study with regard to methods such as LDA+U or 
LDA+DMFT, which supplement a "non-interacting" Kohn-Sham Hamiltonian with a 
Hubbard interaction between a strongly interacting subset of orbitals: using different 
methods for obtaining the noninteracting reference can lead to significant differences, 
and it is by no means clear whether PBE (GGA) or even LDA always provides the best 
starting point for a more sophisticated treatment of correlation effects. Our results also 
emphasize the importance of finding improved ways to account for the double counting 
correction when using different electronic structures as noninteracting reference. 

In order to see how, within the limitations discussed in the preceeding paragraph, 
a TB Hamiltonian of the form ( I6aj) -fl6d j) can reproduce the MLWF band dispersion, we 
consider a modified parameterization using to model the el-el interactions. Since in 
that way the correlation-induced increase of the spin splitting is only partially covered 
by the el-el term ([7]), we correct this by introducing an "empirical" correction to the 
Hund's rule coupling: 

= Jh - Jr"^ - ■ (13) 

Note, that analogously we could choose as the el-el interaction parameter and define 
an appropriate correction to A^. However, since the fundamental band gap in LaMnOa 
is largely controlled by the JT induced splitting between occupied and unoccupied Cg 
bands, and since in a TB model for LaMnOs it seems most desirable to describe the band 
gap correctly, we choose to model the el-el interactions. If the correction AJ-^'' is 
neglected, the local majority spin bands around the band gap are still described quite 
well, even though the splitting with respect to the local minority spin bands will be 
underestimated, which might be acceptable for certain applications. 

Figure [9] also shows the dispersion calculated from such a modified TB model with 
explicit el-el interaction, where the correlation induced change of the spin splitting and 
band gap is described by two interaction parameters, and AJ^\ while Jh, A''^, A-*-, 
and rjx are fixed at their respective PBE values. In addition, the hopping amplitudes are 
set to the values given in table [3l The band dispersions using these sets of parameters 
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(shown as green crosses in figure [H]) again almost perfectly follow the MLWF bands. 
The agreement between the bands calculated within the two parameterizations (Model 
1 and 2) also reflects the transferability of the on-site parameters between the structures 
with and without the GFO distortion. 

4. Summary 

In this paper we have presented a general scheme to parameterize, within a TB 
picture, the band structure of the prototypical JT distorted Cg perovskite LaMnOs 
by means of a suitable downfolding of the ab initio electron dispersion relations onto a 
small set of MLWFs. The tabulated TB parameters should provide an interpretative 
direction for more sophisticated many-body model Hamiltonian investigations of similar 
systems [H ESI EZl [78] . 

By comparing the PBE and beyond-PBE findings we can draw the following 
conclusions: 

(i) Ab initio electronic structure results. We find that all methods consistently 
find a Mott-Hubbard insulating state. GWq provides the best agreement with 
experiments in terms of bandgap value, and both GWq and HSE convey a 
satisfactory description of valence and conduction band spectra. While in the 
PBE+[/ and HSE cases a suitable adjustment of the parameters U and Omix can 
selectively improve the performance with respect to either bandgap or magnetic 
exchange interactions, a universal value that provides all quantities with good 
accuracy cannot be found. Even though the standard value amix = 0.25 in HSE 
seems to provide rather accurate magnetic coupling constants, clearly a smaller Omix 
is necessary to obtain a better Mott-Hubbard gap. While the two different codes 
used in the present study lead only to marginal differences in the Kohn-Sham 
band structure and the corresponding TB parameterization, the relative energies 
of different magnetic configurations depend on subtle details of the used methods, 
which hampers a concise comparison between the different energy functional (it 
should be noted however, that the PAW approach is usually considered superior to 
pure pseudopotential schemes). Within VASP a value for the Hubbard U between 2- 
3 eV leads to similar magnetic coupling along c as HSE, but somewhat stronger FM 
coupling within the ab planes. Despite all its well-known limitations when applied 
to strongly-correlated materials, PBE does not seem to perform too badly (of course 
the fact the we have used the experimental structure helps in that respect, since 
PBE is known to fail in properly reproducing the JT distortion in LaMnOa [60]). 

(ii) MLWFs. Despite the difficulties to fully disentangle the effective Cg bands from 
other bands with similar energies, which are most pronounced within PBE+f/ and 
HSE, the resulting MLWFs and associated ordering (Fig|7]) look rather similar and 
are in good agreement with the precedent plots of Yin[13j. This represents a further 
proof of the quality and reliability of the wannier construction of the Cg {Sz"^ — r^) 
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and — y"^) orbitals. Despite these similarities, the differences in the underlying 
band structures lead to distinct differences in the Hamilt nian matri xelements in 
reciprocal space, and allow for an accurate quantitative analysis of the differences 
between the various approximations for the exchange-correlation kernel. 

(iii) TB parameterization. We have demonstrated that the methods-derived changes in 
the TB parameters due to the different treatment of the el-el exchange-correlation 
kernel in conventional and beyond-PBE approaches can be accounted for using 
two different routes: (a) Model 1 {Htb = i^kin + -^Hund + -f^jx)- In this model 
the TB Hamiltonian does not explicitly incorporate an el-el interaction term. All 
changes in the beyond-PBE band structure with respect to the "noninteracting" 
PBE bands are integrated in the hopping, JT and Hund parameters (in particular 
t^"^, A^, and Jh). (b) Model 2 (^tb = ^kin + ^Hund + ^^jt + ^e-c)- In this second 
type of parameterization we have build in an el-el term in the TB Hamiltonian 
explicitly. The el-el interaction effects are treated by parameterizing the on site 
Hund and JT parameters into a noninteracting (PBE) and interacting (dependent 
on u!^^ and U^^) part. Since we found that ^ in order to achieve a 

correct parameterization it is necessary to fix one ?7w channel {U^) and evaluate 
the changes on the remaining one (Aj|^^). Both, Model 1 and 2, yield excellent 
TB bands, essentially overlapping with the underlying MLWFs ones. 
We note that the different levels of approximation for the non-interacting band 
structure can lead to significant changes in the hopping amplitudes, which cannot 
easily be accounted for by a local double-counting correction. In addition, we 
have also shown that the influence of the beyond-PBE treatment on the model 
parameters of the local Hamiltonian cannot be captured by a simple mean-field 
Hubbard term with only one interaction parameter. For an accurate many-body 
or effective model treatment of LaMnOs and similar materials it thus seems most 
desirable to start from the most realistic single particle band-structure (i.e. not 
necessarily LDA or GGA) and use an appropriate double counting correction. The 
exact form of such a correction term, however, is still unclear at this point. A 
possible alternative solution to correctly treat correlation effects without being 
contaminated by the double-counting problem is the GW-I-DMFT scheme, which 
has recently attracted several research groups and will be most likely available in 
the next future [79118(1]. 

In summary, we have shown that MLWFs can be constructed efficiently not only at 
conventional DFT level (PBE), but also from hybrid functional (HSE) and quasiparticle 
(GWo) wavefunctions, through the creation of an appropriate interface between the 
electronic structure code VASP [l2lll3] and the publicly available WannierQO code [53] . 
Thereby, we have used the well-established PW2WANNIER90 interface as benchmark 
at the PBE and PBE-|-f/ level [39]. Given the booming application of hybrid DFT and 
GWo calculations for a wide variety of materials for which the possibility to describe the 
relevant physics using a minimal basis set is important (these include, e. g., Fe-based 
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superconductors [SI], cuprates [Sg and multiferroics [531 El]), the VASP2WANNIER90 
interface which allows to construct MLWFs directly from the widely-used VASP code, 
will provide a valuable tool for future research. From the practical point of view, we 
have demonstrated that MLWFs can be efficiently used to accurately interpolate the 
HSE and GWq band structure from the coarse uniform k-points mesh to the desirable 
(and dense) symmetry lines, thereby remedying the fundamental practical limitation of 
HSE and GWq scheme in computing energy eigenvalues for selected k-points [521 185]. We 
expect that our study will serve as a reference for future studies involving MLWFs-based 
downfolding procedure. 
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Appendix: Tight binding model parameterization 

In this section, we supply exact definitions for all parameters included in our TB models 
and describe how they are obtained from the MLWF Hamiltonian matrix. 

In the following, we express the hopping amplitudes to-,a(K+Ai?)6(fi) as 2 x 2 matrices 
with respect to the orbital indices: t^^a{R+AR)b{R) — ^ 1'''^' (Ai?). There are 4 different 
hoppings parameters which we consider here, the spin-dependent nearest-neighbor 
hopping t***, the second- nearest neighbor hopping t^^ and the second- nearest neighbor 
hopping along the x, y, z axes t"^^. The nearest-neighbor hopping matrix is expressed as 

f^'(±i) = -it^^'(l + r^), (1.1a) 

t;'{±x) = -\t''{2 ■ 1 - 73 • r^' - r") , {Lib) 

(and analogously for tf^'(±y), see [31]). The t^^ which determines hopping amplitudes 
within FM ordered planes is via fll.l^l) (from the Hamiltonian matrix elements for the 
JT{Q^) distorted structure) given as 

f^=(X-|/^^2)^ (1-2) 

In the GEO distorted structure, the nearest-neighbor hopping amplitudes are reduced 
by a factor {1 — rj^), where the coefficient rj^ is calculated as 

(13) 

'* t-[JT(Q^)] ■ ^ ^ 

The hopping parameters are denoted by the corresponding crystal structure (in square 
brackets), for which they are calculated. The t'^'^ parameter, which determines the 
hopping amplitude between A-AFM ordered planes, is taken as an average of t^^ and 
xhe JT-induced splitting of the nondiagonal elements of the hopping matrix within 
the xy plane is incorporated as an additional contribution to the in-plane hopping 

Ar(±i=) = AQ^(i-I^), (1.4) 

(and analogously for At{±y)). The A parameter is determined for the JT{Q^) distorted 
structure as 

A = ^ {Uhu - hl^f + \{hl, - hl,f) , (1.5) 

and is also reduced by the (1 — rjl) factor in the GFO distorted structure. Both, the 
second nearest-neighbor hopping t^^ and the second- nearest neighbor hopping along the 
X, y, z axes are determined (from the JT(Q^) distorted structure) as 

t^' = -\miy + {h-^)'] , (1.6a) 

= -imiY + {h\iY] . (1.65) 

The matrices related to the t^^ hopping parameter are then given by (see e.g. [50] ) 

t{±x ±z) = -t^-^(-l + 73 . r^- - i") , (1.7a) 
t(±a;±j/) = -r^(-l +2.r"), (1.76) 
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(and analogously for t{±y ± i))- The matrices related to the hopping parameter 
have the same form as fll.lal) and fll.lfel) . In the GFO distorted structure, the matrix 
elements are also reduced by (1 — r^f). 

The Hund's rule coupling strength is determined using (l66j) from the orbitally 
averaged spin splitting of the Hamiltonian on-site diagonal matrix elements (for the 
Pbnm structure) 

JB = l[{h'u + hl,Y-{h'n + hl2y] . (1.8) 

The (generally spin-dependent) JT coupling parameter A* is determined (from the 
JT{Q^) distorted structure) as 

where the JT induced eigenvalue splitting Ae of the Cg subspace 2x2 on-site matrix is 
calculated as 



Ae 



{h',,-hi,Y + {2hi,YY^\ (1.10) 



The JT coupling A'^ is effectively reduced due to the GFO distortion mode (see [31] for 
more details) via the rjx parameter calculated as 

„ Ae^Pbum] |Q[JT(Q-)]| 

A£t[JT(g-)] \Q[Pbnm]\ ' ^ ' ' 

where \Q\ = ^/WVTW?- 
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